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Abstract. A modified perturbation theory in the strength of the nonlinear 
term is used to solve the Nonlinear Schrodinger Equation with a random po- 
tential. It is demonstrated that in some cases it is more efficient than other 
methods. Moreover we obtain error estimates. This approach can be useful 
for the solution of other nonlinear differential equations of physical relevance. 



1. Introduction 

We consider the problem of dynamical localization of waves in a Nonlinear 
Schrodinger Equation (NLSE) [1] with a random potential term on a lattice: 

(1.1) id t tp = -J [i> (x+l) + i> (x - 1)] + e x ip + p |V| 2 V, 

where ip = ip (x, t) , x € Z; and {e x } is a collection of i.i.d. random variables 
uniformly distributed in the interval [— -y-, • 

The NLSE was derived for a variety of physical systems under some approxima- 
tions. It was derived in classical optics where tp is the electric field by expanding 
the index of refraction in powers of the electric field keeping only the leading non- 
linear term [2]. For Bose-Einstein Condensates (BEC), the NLSE is a mean field 
approximation where the term proportional to the density /3\ip\ 2 approximates the 
interaction between the atoms. In this field the NLSE is known as the Gross- 
Pitaevskii Equation (GPE) [3, 4, 5, 6, 7, 8]. Recently, it was rigorously established, 
for a large variety of interactions and of physical conditions, that the NLSE (or the 
GPE) is exact in the thermodynamic limit [9, 10]. Generalized mean field theo- 
ries, in which several mean-fields are used, were recently developed [11, 12]. In the 
absence of randomness (1.1) is completely integrable. For repulsive nonlinearity 
([3 > 0) an initially localized wavepacket spreads, while for attractive nonlinearity 
(/3 < 0) solitons are found typically [1]. 

For (3 = this equation reduces to the Anderson model [13], 

(1.2) idtil) = - J [ip (x + 1) + ip {x - 1)] + e x ip. 

It is well known that in ID in the presence of a random potential with probability 
one all the states are exponentially localized [13, 14, 15, 16]. Consequently, diffusion 
is suppressed and in particular a wavepacket that is initially localized will not spread 
to infinity. This has been very recently extended to the many-body particle system 
[17, 18, 19]. This is the phenomenon of Anderson localization. In 2D it is known 
heuristically from the scaling theory of localization [20, 15] that all the states are 
localized, while in higher dimensions there is a mobility edge that separates localized 
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and extended states. This problem is relevant for experiments in nonlinear optics, 
for example disordered photonic lattices [21], where Anderson localization was found 
in presence of nonlinear effects as well as experiments on BECs in disordered optical 
lattices [22, 23, 24, 25, 26, 27, 28, 29, 30]. The interplay between disorder and 
nonlinear effects leads to new interesting physics [28, 29, 31, 32, 33, 34]. In spite of 
the extensive research, many fundamental problems are still open, and in particular, 
it is not clear whether in one dimension (ID) Anderson localization can survive the 
effects of nonlinearities (see however [35, 36]). 

A natural question is whether a wave packet that is initially localized in space 
will indefinitely spread for dynamics controlled by (1.1). A simple argument indi- 
cates that spreading will be suppressed by randomness. If unlimited spreading takes 
place the amplitude of the wave function will decay since the L 2 norm is conserved. 
Consequently, the nonlinear term will become negligible and Anderson localization 
will take place as a result of the randomness. Contrary to this intuition, based on 
the smallness of the nonlinear term resulting from the spread of the wave function, 
it is claimed that for the kicked-rotor a nonlinear term leads to derealization if it 
is strong enough [37]. It is also argued that the same mechanism results in der- 
ealization for the model (1.1) with sufficiently large (3, while, for weak nonlinearity, 
localization takes place [37, 38]. Therefore, it is predicted in that work that there is 
a critical value of (3 that separates the occurrence of localized and extended states. 
However, if one applies the arguments of [37, 38] to a variant of (1.1), results that 
contradict numerical solutions are found [39, 40]. Recently, it was rigorously shown 
that the initial wavepacket cannot spread so that its amplitude vanishes at infi- 
nite time, at least for large enough (3 [41]. It does not contradict spreading of a 
fraction of the wavefunction. Indeed, subdiffusion was found in numerical experi- 
ments [37, 41, 42]. In different works [42, 43, 44] sub-diffusion was reported for all 
values of (3. It was also argued that nonlinearity may enhance discrete breathers 
[33, 34]. In conclusion, it is not clear what is the long time behavior of a wave 
packet that is initially localized, if both nonlinearity and disorder are present. The 
major difficulty in numerical resolution of this question is integration of (1.1) to 
large time. Most researchers who run numerical simulation use a split-step method 
for integration, however it is impossible to achieve convergence for large times, and 
therefore some heuristic arguments assuming that the numerical errors do not af- 
fect the results qualitatively, are utilized [37, 43]. However it is unclear whether 
those arguments apply to (1.1). The motivation of the current work is to propose 
a numerical scheme based on a modified perturbation theory developed in [45, 46] 
which will allow integration of (1.1) and similar equations up to large times and 
with some control of the error based on the form of the remainder term obtained 
in [46]. 

The advantage of the perturbative method is that it provides an estimate of the 
error while there is no such estimate in the split-step method. Moreover the error 
in the split-step method is expected to proliferate as a result of the nonlinearity. 

In Section 2 we briefly review the perturbation theory developed in [46] . In Sec- 
tion 3 we explain the numerical scheme used to compute the different orders in the 
perturbation theory. In Section 4 we show how the error of the perturbation theory 
could be controlled. In Section 5 we present the comparison between the perturba- 
tion theory and an exact integration. The results are summarized in Section 6 and 
open problems are listed there. 
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In summary, this work demonstrates a numerical implementation of the pertur- 
bation theory for (1.1) in powers of (3 which was described in [46], and evaluates 
its possible use. 

2. The perturbation theory 

Our goal is to analyze the nonlinear Schrodinger equation (1.1) that could be 
written in the form 

(2.1) id t iP = Hoijj + Pm 2 iP 
where H is the Anderson Hamiltonian, 

(2.2) H tp (x) = - J [tp (x + 1) + ip {x - 1)] + e x ip (x) . 

The wavefunction can be expanded using the eigenstates, u m (x), and eigenvalues, 
E m , of H as 

(2.3) V (x, t)=J2 c ™ (*) e ~ tEmtu m (x) ■ 

m 

For the nonlinear equation the dependence of the expansion coefficients, c n (t) , is 
found by inserting this expansion into (2.1), resulting in 

(2.4) id t c n = P J2 Vr im2m3 C i c m2 c ro3 e i ( E «+^ 1 -^-^3)* 

771 1 ,7712 ,7713 

where V™ 1 " 12 ™ 3 is an overlap sum 

(2.5) C" 8 " 13 = E u » ( x ) u ^ ^ u — u -3 (x) ■ 

X 

Our objective is to develop a perturbation expansion of the c m (t) in powers of (3 
and to calculate them order by order in (3. The required expansion is 

(2.6) c n (t) = 4°) + /3c« + (3 2 cW + ■■■+ ^ + Q n , 

where the expansion is till order N and Q n is the remainder term (note here Q n 
differs from one defined in [46]). We will assume the initial condition 

(2.7) c n (t = 0) = S n0 . 

The equations for the two leading orders are presented in what follows. The leading 
order is 

(2.8) = 6 n0 . 
And the first order is 

(2.9) ^=vr{ - ; F 

We notice that divergence of this expansion for any value of (3 may result from three 
major problems: the secular terms problem, the entropy problem (i.e., factorial 
proliferation of terms), and the small denominators problem. In this paper we will 
not discuss the entropy problem and the problem of small denominators, since it 
was done in detail in [46]. We first show how to derive the equations for c„ (t) 
where the secular terms are eliminated. To achieve this we replace the ansatz (2.3) 

by 
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(2.10) V (x, t) = Y^Cn (t) e- lE '^u n (x) 

n 

where 

(2.11) K = + + (3 2 E^ + ■■■ 

and E^ are the eigenvalues of H . We will dub E' n the renormalized energies. The 
new equation for the c n is given by 

(2-12) 

id t c n = (4°) - <) c n + P Yl ^ miTO2ro3 ^ 1 c m2 c m3 e<< +B -- B --^3)*. 

Inserting expansions (2.6) and (2.11) into (2.12) and comparing the powers of (3 
without expanding the exponent in /?, produces the following equation for the k — th 
order 

(2-13) 



fc-i 



s=0 

'k-1 k-l-r k-l-r- 



EE E 

.r=0 s=o l=o 



i(E' +E' —El -E', )t 
\ n mi m 2 »3 ; 



Note that the exponent is of order O (1) in /?, and therefore we may choose not to 
expand it in powers of (3. However, it results in an expansion where both Em and 

(k) (I) (k) 

Cn depend on f3. For the expansion (2.6) to be valid, both E m and c„ should 
be O (1) in (3, this is satisfied, since the RHS of (2.13) contains only c£ r ' ) such that 
r < k. Namely, this equation gives each order in terms of the lower ones, with the 
initial condition of ci°^ (t) = <5„o ■ Solution of k equations (2.13) gives the solution 
of the differential equation (2.12) to order k, while the higher order terms which 
are obtained from this equation are meaningless (see Appendix for the reasoning). 
Since, the exponent in (2.13) is of order O (1) in (3 we can select its argument to 
be of any order in (3. However, for the removal of the secular terms, as will be 
explained bellow, it is instructive to set the order of the argument to be k — 1, 
as the higher orders were not calculated at this stage. Secular terms are created 
when there are time independent terms in the RHS of the equation above. We 
eliminate those terms by using the first two terms in the first summation on the 
RHS. We make use of the fact that ci°^ = S n o and Cn ^ can be easily determined 
(see (2.16,2.15)), and used to calculate E^2 and -Ew that eliminate the secular 



.(*) 

-•71 ■ 

(2.14) 4 fc )cl°) + E^^ = E^S n0 + E^ (1 - M 



terms in the equation for c„ , that is 



K-K' 



where only the time-independent part of Cn^ was used. In other words, we choose 
E^ and E^^ so that the time-independent terms on the RHS of (2.13) are 

(k) (k—1) 

eliminated. E will eliminate all secular terms with n = 0, and En will 
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eliminate all secular terms with n ^ 0. In the following, we will demonstrate this 
procedure for the first order. 

In the first order of the expansion in f3 we obtain 

(2.15) idttf = -EMcM + J2 y ri -=-3 c Ho) c (o) c (o) e <K + ^ 1 -^ 2 -^ 3 )' 

771177127713 

= -E^6 n0 + V° 00 e< K - E ' o)t - 
For n = the equation produces a secular term 

(2.16) idJV = -E™ + V a oao 



(2.17) = V 000 



Setting 
(2.17) 

eliminates this secular term and gives 

(2.18) c 1} = 0. 

For n 7^ there are no secular terms in this order, therefore finally 

'l_ e «-^)t\ 



(2.19) c« = (1 - S n0 ) V 7 



0(11] 



E' n - E' 



where to this order E' n = E n and E' = Eq. 

The higher order terms in the perturbation theory are given by recursive rela- 
tions and due to the large number of terms which are involved will be calculated 
numerically. 

3. The numerical method 

In order to compute the various orders in the perturbation theory we use equation 
(2.13), which is a recursive equation of the orders. To compute order k we have 
to compute all cl and E„ for I < k — 1. The numerical calculation is done in 
two stages: at the first stage a symbolic calculation of the expressions of all the 
Cn and En is performed, this has a complexity of O (e 2k ) , which is due to the 

increasing number of terms in each expression for e'n (see [46]). This stage does not 
depend neither on the realization nor the nonlinearity strength, (3. In the second 
stage realizations and (3 are chosen and c« and E® are calculated. This stage 
has a complexity of O (e 2fc • L fe ) , where L is the dimension of the lattice. The 
computation of this stage could be fully parallelized. 

When calculating E n ^ we encounter self-consistent equations of the type 

(3-1) E n = f n ({<„}) , 

where / is some function, for example for the second order 



(3.2) E' n = 4°) + /3V/ n " 00 (2 - S n0 ) - 3(3 2 S n0 



Higher order equations are required in general. We solve those equations numeri- 
cally by reinserting the LHS into the RHS, until a desired convergence is achieved. 
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The first iteration is done by setting (3 = at the RHS. Basically, at each iteration 
an order of (3 is gained in the accuracy of the solution and since we need to know 
E n only to a desired order N (see Appendix) , only a small number of iterations is 

needed. The c$ are represented as vectors with elements fc^i, c^;L 2 , ■ ■ -J identi- 
fied by frequencies such that terms with same frequencies are grouped together (by 
summing their amplitudes), namely, Cn^ k = J2j c nL fc j e ~ tIJt 'i where Uk is a shared 
frequency. Due to the fact that most of the amplitudes are negligible, after grouping 
a thresholding step is done and terms which are smaller than 10 -6 are eliminated. 
The error introduced by the tresholding can be easily controlled, since we know 
how many frequencies were left out. By having the vector of frequencies and their 
corresponding amplitudes we can calculate the perturbative solution at any time. 
Even after grouping and thresholding the number of frequencies is growing rapidly 
with the order of the expansion. 



4. The remainder of the expansion 



In order to control the solution we have to control, Q n , the remainder of the 
expansion (2.6) that can be written in the form 

(4.1) c„ (t) = c„ + Q n , 

with 

N 



(4.2) 

It is useful to define 
(4.3) 

and 
(4.4) 



5> 

(=0 



(0 



-iE„t 



Qn Qn& 



-iEt 



Substituting (4.1) in (2.12) leads to the following equation for the remainder which 
is expressed in terms of (4.2), (4.3) and (4.4), 

(4.5) ld t Q n = W n (t) + Mnm (t) Qm + £ M nm (t) Q* m + F (q 



where 
(4.6) 



W n (t) = (^ 0) - <) 5 n e- l <* - i (d t Cn) e 



iEt 



2 _ 



+ (3 ^ U n (x) i> (x) 4> (x) 

X 

is the inhomogeneous term, 

2 

(4.7) M nm (t) = E^S nm + 2pJ2 u ^( x ) $( x ) u m{x) 

X 

and 



(4.5 



Mr, 



(t) = ^ u n (x) (j; (x)) u m (a 
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determine the linear terms, while the nonlinear term is, 



(4.9) F(Q) = p^unix)^* {x)\Y j Q m u m {x)\ 

x \ m / 



[3 ^2 u n (x) 



The linear part of (4.5) is given by 



(4.10) 



ld t Ql n = W n (t) + Mnm (*) Q l ™ + Mnm (*) (Q 



Using a bootstrap argument, which utilizes the continuity of (4.5) and smallness 
of the linear part of (4.10) one can show [46] that until some time t* , the dynamics 
of (4.5) is governed by the linear part and the remainder is bounded by, 



(4.11) 



\Q n (t)\ < A-t-e-fW 



where 7 is the inverse localization length. Therefore to estimate the remainder we 
can integrate (4.10) instead of (4.5) at least up to i*. It is useful to integrate up to 
some large time, t <C t*, and then to extrapolate using the linear bound (4.11) up 
to t*. In the next section it will be proposed how to determine t* in practice. 



5. Results 

In this section it will be demonstrated, how the numerical scheme for calculations 
in the framework of the perturbation theory, is implemented in practice. Some 
results will be compared with an exact numerical solution of the original equation 
(1.1). 

For this purpose we have calculated numerically all the c4 and E n for I < 4 
for a certain realization of the random potential. To compare perturbation theory 
results to the exact results, we compute their Fourier transform for different orders 
of expansion. On Fig. 5.1 we see the Fourier transform of Co, c\ and eg (see (4.2)) 
compared to the Fourier transform of an exact (numerical) solution, c n , calculated 
using a split-step method. We notice a reasonable agreement of the perturbation 
theory with an exact solution for cq, c\ and a disagreement for eg . By plotting 
Q l ™ (t) for all n with the same scale (on the same axis) in Fig. 5.2, we see that 
there are modes (on Fig. 5.2 there are two of them) which contribute to most of 
the discrepancy in the perturbation theory calculation, since if Q 1 " 1 (t) is large the 
bound on Q n (t) is also large. 
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0.07 




100 200 300 400 500 600 700 800 900 1000 

t 



Figure 5.2: Q l ™ as a function of time in the 4th order in (3 for all n's of the lattice 
(total 128 lines). The two lines that are far above the rest (which are barely visible) 
correspond to the resonant modes, n = 4,9. The parameters are: (3 = 0.0774, 
W = 4, J = 1. 



We will call those modes resonant modes. In Fig. 5.3 we compare a norm of Q 
calculated with the resonant modes, 

(5.i) iiqii 2 = (£iQ m | 2 

and without them ||Q'|| 2 . 




3 100 200 300 400 500 600 700 800 900 1000 



t 

Figure 5.3: ||<5|| 2 (dashed blue) and ||Q'|| 2 (solid green) as a function of time. The 
straight line is a linear fit to ||Q|[ 2 . The parameters are: 4th order, /? = 0.01, 
W = 4, J = 1. The linear fit is: ||Q|| 2 (t) = 2.3 x 10~ 9 •H5.7x 10~ 6 . 
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It is found that indeed the ||Q|[ 2 grows linearly with time as expected from 
(4.11). Actually, this is the way a bound on the resonant terms is expected to 
behave [46]. It is evident that by excluding the resonant modes the discrepancy 
between the perturbation theory and exact results is much lower. The reason why 
the perturbation theory fails for the resonant modes is that they correspond to a 
quasi-degeneracy, namely, when E' n E' , and the overlap | V® 00 1 is not sufficiently 
small. A natural way to quantify the resonance condition is to use (compare to 
Eq.(5) of [43]), 



(5.2) 



V. 







K - K 



The resonant modes produce substantially higher values of R n 1 compared to any 
other modes as demonstrated in Fig. 5.4. 




48 64 



Figure 5.4: R n x as a function of n 



One way to deal with this problem is by using a degenerate perturbation theory. 
However it is an open question how to implement it for a nonlinear problem, that 
should be left for further studies. 

At the end of the last section we have claimed that for some time, i», the linear 
part of (4.5) dominates over the nonlinear part given that the linear part is suffi- 
ciently small. In Fig. 5.5 we present a comparison between the solution of the linear 
equation (4.10) and Q ex , the solution of (4.5). It is clear that until ||Qim|| 2 ~ 0-1 
both the solutions are very close. 
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Figure 5.5: HQe^^ (lines) and HQ^n H2 (dotes) as a function of time, for various 
values of (3 (see legend). The perturbation expansion is up to fourth order in 



We use this value to define 



(5.3) 



| Qlin (^*) | 



0.1. 



For small nonlinearity strength, (3, t* is very large and therefore the integration 
of (4.10) to is very time consuming. We therefore use the bound (4.11) to 
extrapolate linearly from the time interval where (4.10) is solved to i». Practically, 
we have calculated the linear behavior of Qu n like it was done in Fig. 5.3 than we 



found i* from (5.3). In Fig. 5.6 we plot log 10 i* as a function of (3 



for different 
of 



orders, while in Fig. 5.7 we compare the result found in 4-th order with \Q 1 " 1 
(5.3) is replaced by Q l ™ (where the resonant terms are removed) 
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Figure 5.6: log 10 t* as a function of (3 1 for different orders. 4th order (solid blue), 
3rd order (dashed green) and 2nd order (dotted red). The parameters are: W = 4, 
J = 1. 




^0 20 30 40 50 60 70 80 90 100 

Figure 5.7: log 10 £* as a function of /3 _1 for the 4-th order.. Solid blue line is when 



i* is calculated based on \Q\\ n | and dashed green line is when 
The parameters are: W = 4, J = 1. 



Q'„ 



is used instead. 
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A systematic improvement with the order of the perturbation theory is found. 
We notice that even with moderate nonlinearity strengths, namely, (3 < 0.08, one 
can achieve a good approximation of the solution up to very large times. In this way 
the perturbation theory combined with the solution of the linear equation (4.10) 
and the criterion (5.3) may be used to obtain the solution of the original equation 
(1.1) up to t < t*. For small (3 the time t, is very long as is clear from Fig. 5.6. 
When one considers the smallness of (3 one should consider actually the smallness 
of /3e 2 due to the exponential proliferation of the number of terms (see Eq. 4.6 of 
[46]). 

Assuming that the Q' n give a good approximaion of the Q n for t < £* , we can 
use Fig. 5.7 to conclude that, for example, for j3 = 0.1 our solution is meaningful 
up to time t = 10 4,3 and therefore can be whithin reach of works like [43, 44]. 
The assumption obove corresponds to the observation that correcting C4 and c g by 
Q4 and Qg is not changed much over time, and therefore all the Q n stay small 
for a very long time. Note, that since the largest localization length for W — 4 
is approximately 6, positions 4 and 9 are within one localization length from the 
initial data. Therefore their removal does not affect the behavior for large n that 
is relevant for the asymptotics. 

Additional benifit of this implemamtnation is the computational speed. For 
example, even for the 3rd order of the pertrurbation theory one can compute the 
perturbative solution up to t < 10 4 in 10 minutes while using the same computer 
an exact convergent integration (split-step method) takes 28 hours, which is a two 
orders of magnitude speed-up. 

6. Summary and Discussion 

In this paper we have demonstrated how a perturbation theory can be numeri- 
cally implemented for the solution of the NLSE with a random potential (1.1). In 
the perturbative method the computer is used to implement the symbolic recur- 
sive calculation of the various terms and also for their numerical evaluation. This 
method allows to estimate the errors since the remainder term is bounded. It was 
demonstrated (in the end of the previous section) how evaluate the time of validity 
of the perturbation theory. This approach has also a great advantage in the speed of 
the calculation. We believe that the method can be generalized to other nonlinear 
differential equations, for example, the Fermi-Ulam-Pasta problem. 

In order to make the perturbation theory of a much greater value one has to solve 
the problem of resonant terms (when i?" 1 of (5.2) is large). In other words the 
degenerate perturbation theory should be extended to nonlinear equations. Also 
the asymptotic nature of the perturbation theory is not known, it is not clear when 
it is convergent and when only asymptotic. 
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research fund. 
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Appendix 

In the calculations of the paper we used the expansion (2.10) of the wavefunction 
and the expansion coefficients c„ (t) are calculated from (2.6) and (2.13). The ci (t) 
depend on j3 since the E n in (2.13) depend on it. Therefore (2.6) is not a Taylor 
series in (3. In this Appendix it is shown that this expansion is equivalent to the 
Taylor series, 

oo 

i=o 

Contrary to the cl the c$ are independent of /3. The expansion of (2.12) in powers 
of p produces the following equation for the r-th order 

(6.1) 

1 fd h E { t h) \ 



••*4 r) = -EEm-^5— j c 

1=1 h=0 \ ' / f3 =0 



(r-l)_ 



V V ymim2m 3 ,.(!i)=(i 2 ),(! 3 ) ( J_ SL!_j{K+E' mi -E' -E' )t 
ZL/ ™ till u m 2 °m 3 I m( 4 g 



E,i.=''-l{m i } 



We considered a different expansion (2.13), where one does not expand the expo- 
nent, but still compares the implicit powers of (3 from both sides of the equation. 
The equation for the (r — k) order in this expansion is (that is just (2.13)) 

(6.2) 

1=1 

, V V V «iirn 2m3 *(h) (l 2 ) <l 3 ) i(E' n +E' mi ~E' mo -E' )t 

' / j / j v n u mi °rri2 u m3 Q 

{mi} Ii+l2+l3=r — k — l 

It can be shown that each term in this expansion is uniformly bounded in time, 
with a proper selection of E' n , that is the secular terms are removed. We will show 
that this expansion is equivalent to a Taylor series. 

Theorem 1. For any t, 

oo 

(6.3) c n =Y,c { n ) P k , 

k=0 

where c4 are defined using (6.2). 

Proof. Since c n = J2T=o (3 k is the expansion of c n in powers of /3, and due to 
the uniqueness of this expansion, the theorem is true if and only if 

k=0 V ^ / 0=0 k=0 \ / /3=0 

where at the last equality we have used the fact that d g 'j 3 " = 0, for any r > 1. We 
proceed to prove this equality by induction. Suppose that this equality is true for 
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all the orders till order r. We apply the linear operator J2k=o to both sides 

of (6.2). For the first part on the RHS we get 



r— 1 -, r—k ofc 

f^k\4-^ d0 k V " 



(r-fc-0 



fe=0 1=1 



^ fc! ^ h\(k-h)\ df3 h 5/3( fe -'i) 



(6.5) 



r-lr-fe fe 

EEE 

fe=0 Z=l d=0 



1 d' 1 ^ 



(Jfc-li)! <9/3( fc -'i) 



where at the first equality we have used the Leibniz generalized product rule. We 
now exchange variables such that 



(6.6) 

or in a matrix notation 

(6.7) 

where 



h = h 
h = k - 1\ 
z = l + h 



h 
h 

z 



h 

= A 1 \ I 



1 
A x = | -1 1 
1 1 

The region of summation is bounded by the planes 



(6.8) 



or in matrix notation 



< k < r - 1 

1 < l<r-k 
< h < k 



(6.9) 



where 



Bx 



h 
I 
k 



( r-1 \ 
r 



-1 

V o J 



Bx = 



{ 1 \ 

1 1 

1 -1 
-10 

0-10 

V o o -l J 
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Changing the variables by applying the transformation matrix produces 

/ r-n 

r 



-1 

V o j 



(6.10) 




( 


1 


1 





\ 







1 


1 









-1 









-1 












1 





-1 




\ 


-1 


-1 





) 




< 



or the following inequalities 
(6.11) 



h+h < r-1 

h + z < r 

< h 

< h 

h < z-1 

< h+h- 
Removing redundant inequalities gives 

(6.12) < h < z-1 

1 < z <r 
< h<r~z. 

Therefore in the new variables the sum (6.5) takes the form 



(6.13) 



Since 



1 d ll E, 



hi 



dp* 



(6.14) 



r z — 1 r — z 

EEE 

z=i i 1= o ; 2 =o 
for any z > l,we can write 
1 d h Ei z ~ hY 



1 d 



hi 



d(3 l 



z r — z 



EEE 

z=l i 1= 0Z 2 =0 



1 d 1 '^ 



hi dp* 



hi 



d(3 l 



Taking (3 = and using the assumption of the induction for orders lower than r we 
have 

' 1 d h E [ n z - ll)S 



(6.15) 



EE 

z=l l 1= o 



hi dp* 

5 W 



which is the first expression for idtc n ■ The proof for the second term is similar, 

,r-l 1 d k 
<k=0 k\ d/3 k 



operating with X)fe=o F^F on the second term in (6.2) gives 



(6.16) 



r — 1 r — k—1 r — k—l — li 



£v~ m3 E E E 

{mi} k=0 h=0 Z 2 =0 



1 d k 



r .*(h) r (l2)Jr~k-l-h-h) AE' n +E' -E' -E' )t 
L mi 7712 ^ 



Using the Leibniz generalized product rule, which states 

d k . . v-^ kl d^x 1 d S2 x 2 d S3 x 3 d^x 4 

{x lX2 x 3 x 4 ) = ^ 



( 6 - 17 ) ~^k { x l x 2XsX A ) 



si+s 2 +s 3 +Si=k 



si!s 2 !s 3 !s 4 ! <9/3 Sl df3 s * df3 s => d(3 s * 
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we have the sum 



r-lr-fc-lr-fc-l-Ji k k-si k-si-a-, / oai / i as 2 Jh) 

ym 1 m 2 m3 ^ 

{mi} fe=0 h=0 l 2 =0 si=0s 2 =0 s 3 =0 

(6.18) 



si! 9/3 s i / ls 2 ! 9/3 

1 d' !3 C m '~ fc_1_ ' 1_i2 ^ \ / 1 g(*-ai-S2-«a) 



s 3 ! 9/3 s 3 / - sj - s 2 - s 3 )! 9 ( 5(fe-*i-s2-s 3 ) 



«+< I -B^ 3 -< s )* 



Following the first part of this proof we exchange the variable to 



(6.19) 



Z\ 

Z3 
Si 



h + si 

^2 + S 2 

fc — Si - s 2 
Si 



S3 



or using a transformation matrix 



(6.20) 



( Z1 ) 




( k \ 


Z2 




h 


Z3 


= A 2 


h 


Si 




Sl 


S2 






\ S 3 j 




I S 3 j 



where 



A 2 = 



( 


1 





1 
















1 





1 







1 








-1 


-1 


-1 













1 






















1 







I o 














1 


) 



The region of summation is bounded by the following hyperplanes 



(6.21) 






< 


k < r - 


1 







< 


h<r- 


k- 


1 





< 


h<r- 


k- 


1-Zi 





< 


Si < k 









< 


s 2 < k - 


- Sl 







< 


s 3 < k - 


- Sl 


- s 2 
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which could be represented in a matrix notation as 



(6.22) 



B 2 



( k \ 

h 
h 
si 

V ^3 J 



< 



( r-1 \ 
r - 1 
r - 1 





















where 



B 2 = 
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1 
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-1 
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1 





-1 
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1 
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-1 




















-1 




















-1 























Using the transformation matrix to change the variables results in the following 
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or in the following inequalities 
(6.24) 



Z 3 " 


h Si - 


h s 2 - 


f S3 
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r — 


1 




z\ - 


hz 3 " 


h S 2 - 


f S3 
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r — 


1 




z\ - 


V Z-2 - 


h2 3 - 


f S 3 


< 


r — 


1 













< 


23 - 


hs 2 - 


-S3 











< 


z 3 - 


hs 3 













< 


Z 3 















< 


z-i - 


V si - 


- s 2 








Si 


< 


Z\ 












S2 


< 


?2 















< 


Si 







S3 



Removing redundant inequalities gives 

(6.25) < si < 2i 

< s 2 < z 2 

< s 3 <(r-l)-zi-z 2 - z 3 

< zi < (r - 1) - z 3 

< z 2 < (r - 1) - z x - z 3 

< z 3 < r - 1 

which is equivalent to the sum 

r— 1 7 — l-z 3 r — l-zi— z 3 zi z 2 r-1 — Z1-Z2-Z3 / -, a Sl *(zi— si)' 

E v mim 2 m 3 \ " \ " \ " \" \" \" [ 1 °' Cm i 

" Z^ I s ! 5«s! 

{m,} z 3 =0 zi=0 z 2 =0 si=Os 2 =0 s 3 =0 \ 

(6.26) 

/ 1 a s2 c^r S2) \ / 1 d^ct' 1 '" 1 '" 2 '^ 33 ^ ( 1 d z * 




s 2 \ 9/3 S2 / \s 3 ! d(3 s * / 



i(E'+E' -E' -E' )t 

\ 71 77i 2 m 3/ 



Putting /3 = and utilizing the assumption of the induction we get 
(6.27) 

V V '\/ m i m 2 m 3 p(h) ( _LJ^XK+E' mi -E' mo -E ' )t 

Z^ Z^ V » C ™2 C ™3 I l l df3 U e 

li+h+h+U=r-l { mi } v ' 7 P=° 

_(r) 

which is exactly the second term in the equation for idtc„ . Since the zero order 
trivially satisfies this theorem, this competes the proof by induction. □ 

One should be able to extend the results of this Appendix to other nonlinear 
equations, for example, where the power of the nonlinearity is different (in the last 
term of (1.1) 2 is replaced by another integer). 
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to 



Figure 5.1: Fourier transform of cq, C\ and eg for different orders of the perturbation 
theory compared to the Fourier transform of an exact solution, c„, (solid line). 
Second order is given by green crosses, third order by red triangles and fourth 
order by blue squares. The parameters are: (3 = 0.0774, t = 1000, W = 4, J = 1. 



